#include "argList.H"
#include "Time.H"
#include "fvMesh.H"
#include "fvc.H"
#include "fvm.H"
#include "timeSelector.H"
#include "labelIOField.H"
#include "labelFieldIOField.H"
#include "labelList.H"
#include "scalarIOField.H"
#include "scalarFieldIOField.H"
#include "vectorIOField.H"
#include "vectorFieldIOField.H"
#include "vectorList.H"
#include "tensorIOField.H"
#include "tensorFieldIOField.H"
#include "tensorList.H"
#include "pointFields.H"
#include "fluidThermo.H"
#include "fluidThermoMomentumTransportModel.H"
#include "fluidThermophysicalTransportModel.H"
#include <iostream>
#include <fstream>
#include <vector>
#include <array>
#include <cmath>
#include <numeric>
#include <limits>
#include <string>
#include <exception>
#include <Eigen/Dense>

/* Constants */
const Foam::scalar s2ns = 1e9, ns2s = 1.0/s2ns;
const Foam::scalar s2us = 1e6, us2s = 1.0/s2us;
const Foam::scalar s2ms = 1e3, ms2s = 1.0/s2ms;
const Foam::scalar m2nm = 1e9, nm2m = 1.0/m2nm;
const Foam::scalar m2um = 1e6, um2m = 1.0/m2um;
const Foam::scalar m2mm = 1e3, mm2m = 1.0/m2mm;
const Foam::scalar m2cm = 1e2, cm2m = 1.0/m2cm;
const Foam::scalar kcal2J = 4186.8;
const Foam::scalar one_atm = 101325.0;           // Unit: Pa
const Foam::scalar one_bar = 100000.0;           // Unit: Pa
const Foam::scalar one_psi = 6894.76;            // Unit: Pa
const Foam::scalar one_mpa = 1e6;                // Unit: Pa
const Foam::scalar G0 = 1.4;                     // Specific heat ratio for ideal gas

// Universial gas constant, Unit: J/mol/K
const Foam::dimensionedScalar R(Foam::dimEnergy/Foam::dimMoles/Foam::dimTemperature, 8.31446261815324);

// Molecular weight, Unit: kg/mol
const Foam::dimensionedScalar MW(Foam::dimMass/Foam::dimMoles, 26e-3);

// Gas constant, Unit: J/kg/K
const Foam::dimensionedScalar Rg(R/MW);

/* Classification of gas-phase cells */
const Foam::scalar cFluid = 1.0;                 // Fluid cell
const Foam::scalar cGhost = 0.0;                 // Ghost cell
const Foam::scalar cSolid = -1.0;                // Solid cell

const Foam::scalar cFresh = 1.0;

/* Flow condition */
const Foam::scalar Pr = 0.71;                    // Prandtl number
const Foam::scalar p0 = 2.07 * one_mpa;
const Foam::scalar T0 = 300.0;

/* Plane param */
const Foam::scalar plane_z = 0.25e-3;            // Vertical position
const Foam::scalar plane_T = 800.0;              // Temperature
const Foam::scalar plane_rb = 100 * mm2m;        // Escape velocity, Unit: m/s

/* Cartesian grid */
const Foam::scalar L = 0.5e-3;                   // Domain length, Unit: m
const Foam::label N = 64;                        // Grid resolution
const Foam::scalar h = L / N;                    // Spacing
const Foam::scalar h_inv = 1.0 / h;
const Foam::scalar xMin = 0.0, xMax = L;         // Range in X-direction
const Foam::scalar yMin = 0.0, yMax = 0.1*L;     // Range in Y-direction (empty)
const Foam::scalar zMin = 0.0, zMax = 2*L;       // Range in Z-direction

/* Connectivity of grid cells (degree 2) */
struct adj2Stencil
{
    Foam::label myC;
    Foam::labelList adjC;
};

inline bool isEqual(double x, double y)
{
    return std::islessequal(std::abs(x-y), 1e-5 * std::abs(x));
}

inline bool isZero(double x)
{
    return std::isless(std::abs(x), 1e-7);
}

inline void xyz2ijk(const Foam::vector &p_, int &i_, int &j_, int &k_)
{
    const double dx = p_.x() - xMin;
    const double dy = p_.y() - yMin;
    const double dz = p_.z() - zMin;
    i_ = static_cast<int>(std::fma(dx, h_inv, 0.5));
    j_ = static_cast<int>(std::fma(dy, h_inv, 0.5));
    k_ = static_cast<int>(std::fma(dz, h_inv, 0.5));
}

void setBdryVal(const Foam::fvMesh &mesh, const Foam::volScalarField &src1, const Foam::volVectorField &src2, Foam::surfaceScalarField &dst)
{
    for (int pI = 0; pI < mesh.boundary().size(); pI++)
    {
        const auto &curPatch = mesh.boundary()[pI];
        const auto &patch_val1 = src1.boundaryField()[pI];
        const auto &patch_val2 = src2.boundaryField()[pI];
        const auto &patch_Sn = mesh.Sf().boundaryField()[pI];
        auto &patch_dst = dst.boundaryFieldRef()[pI];
        for (int fI = 0; fI < curPatch.size(); fI++)
            patch_dst[fI] = patch_val1[fI] * patch_val2[fI] & patch_Sn[fI];
    }
}

void setBdryVal(const Foam::fvMesh &mesh, const Foam::volVectorField &src, Foam::surfaceScalarField &dst)
{
    for (int pI = 0; pI < mesh.boundary().size(); pI++)
    {
        const auto &curPatch = mesh.boundary()[pI];
        const auto &patch_val = src.boundaryField()[pI];
        const auto &patch_Sn = mesh.Sf().boundaryField()[pI];
        auto &patch_dst = dst.boundaryFieldRef()[pI];
        for (int fI = 0; fI < curPatch.size(); fI++)
            patch_dst[fI] = patch_val[fI] & patch_Sn[fI];
    }
}

void diagnose(const Foam::fvMesh &mesh, const Foam::globalMeshData &meshInfo, const Foam::volScalarField &src, const Foam::volScalarField &flag, double &norm1, double &norm2, double &normInf, double &minVal, double &maxVal)
{
    norm1 = norm2 = normInf = 0.0;
    minVal = std::numeric_limits<double>::max();
    maxVal = std::numeric_limits<double>::min();
    for (int i = 0; i < mesh.nCells(); i++)
    {
        if (isZero(flag[i]))
            continue;

        const double cVal = std::abs(src[i]);
        normInf = std::max(normInf, cVal);
        norm1 += cVal;
        norm2 += cVal*cVal;
        minVal = std::min(minVal, src[i]);
        maxVal = std::max(maxVal, src[i]);
    }
    Foam::reduce(normInf, Foam::maxOp<Foam::scalar>());
    Foam::reduce(norm1, Foam::sumOp<Foam::scalar>());
    Foam::reduce(norm2, Foam::sumOp<Foam::scalar>());
    norm1 /= meshInfo.nTotalCells();
    norm2 = std::sqrt(norm2/meshInfo.nTotalCells());
    Foam::reduce(minVal, Foam::minOp<Foam::scalar>());
    Foam::reduce(maxVal, Foam::maxOp<Foam::scalar>());
}

void diagnose(const Foam::fvMesh &mesh, const Foam::globalMeshData &meshInfo, const Foam::volScalarField &src, const Foam::volScalarField &flag, double &norm1, double &norm2, double &normInf)
{
    norm1 = norm2 = normInf = 0.0;
    for (int i = 0; i < mesh.nCells(); i++)
    {
        if (isZero(flag[i]))
            continue;

        const double cVal = std::abs(src[i]);
        normInf = std::max(normInf, cVal);
        norm1 += cVal;
        norm2 += cVal*cVal;
    }
    Foam::reduce(normInf, Foam::maxOp<Foam::scalar>());
    Foam::reduce(norm1, Foam::sumOp<Foam::scalar>());
    Foam::reduce(norm2, Foam::sumOp<Foam::scalar>());
    norm1 /= meshInfo.nTotalCells();
    norm2 = std::sqrt(norm2/meshInfo.nTotalCells());
}

void diagnose(const Foam::fvMesh &mesh, const Foam::volScalarField &src, const Foam::volScalarField &flag, double &minVal, double &maxVal)
{
    minVal = std::numeric_limits<double>::max();
    maxVal = std::numeric_limits<double>::min();
    for (int i = 0; i < mesh.nCells(); i++)
    {
        if (isZero(flag[i]))
            continue;

        minVal = std::min(minVal, src[i]);
        maxVal = std::max(maxVal, src[i]);
    }
    Foam::reduce(minVal, Foam::minOp<Foam::scalar>());
    Foam::reduce(maxVal, Foam::maxOp<Foam::scalar>());
}

inline bool pntInFluid(Foam::scalar c)
{
    return c > 0.0;
}

inline bool pntInSolid(Foam::scalar c)
{
    return c < 0.0;
}

inline bool pntNearSurf(Foam::scalar c)
{
    return std::abs(c) < 0.1 * h;
}

bool cellIsIntersected(const Foam::scalarList &v)
{
    bool flag = pntInSolid(v[0]);
    for (int i = 1; i < v.size(); i++)
    {
        if (pntInSolid(v[i]) != flag)
            return true;
    }
    return false;
}

Foam::scalar cellNum(const Foam::scalar c, const Foam::scalarList &v)
{
    if (pntInFluid(c))
        return cFluid;
    else
    {
        if (cellIsIntersected(v))
            return cGhost;
        else
        {
            bool awayFromSruf = true;
            for (int i = 0; i < v.size(); i++)
            {
                if (pntNearSurf(v[i]))
                {
                    awayFromSruf = false;
                    break;
                }
            }

            if (awayFromSruf)
                return cSolid;
            else
                return cGhost;
        }
    }
}

void identifyIBCell(const Foam::fvMesh &mesh, const Foam::volScalarField &phi, const Foam::pointScalarField &phi_pnt, Foam::volScalarField &marker)
{
    int nSolid = 0, nGhost = 0, nFluid = 0;
    for (int i = 0; i < mesh.nCells(); i++)
    {
        const auto &incNode = mesh.cellPoints()[i];
        Foam::scalarList val_n(incNode.size());
        for (int j = 0; j < incNode.size(); j++)
            val_n[j] = phi_pnt[incNode[j]];

        marker[i] = cellNum(phi[i], val_n);

        if (isEqual(marker[i], cSolid))
            ++nSolid;
        else if (isEqual(marker[i], cFluid))
            ++nFluid;
        else
            ++nGhost;
    }

    // Ensure all fluid cells not adjacent to solid cells directly (through face) 
    for (int i = 0; i < mesh.nCells(); i++)
    {
        if (isEqual(marker[i], cFluid))
        {
            const auto &cC = mesh.cellCells()[i];
            for (int j = 0; j < cC.size(); j++)
            {
                const auto adjCI = cC[j];
                if (isEqual(marker[adjCI], cSolid))
                {
                    marker[adjCI] = cGhost;
                    --nSolid;
                    ++nGhost;
                }
            }
        }
    }

    // Ensure all solid cells not adjacent to fluid cells directly (through face) 
    for (int i = 0; i < mesh.nCells(); i++)
    {
        if (isEqual(marker[i], cSolid))
        {
            const auto &cC = mesh.cellCells()[i];
            for (int j = 0; j < cC.size(); j++)
            {
                const auto adjCI = cC[j];
                if (isEqual(marker[adjCI], cFluid))
                {
                    marker[adjCI] = cGhost;
                    --nFluid;
                    ++nGhost;
                }
            }
        }
    }

    if (nSolid > 0)
    {
        // Remove redundent ghost cells on fluid side (through face connection)
        for (int i = 0; i < mesh.nCells(); i++)
        {
            if (isEqual(marker[i], cGhost))
            {
                const auto &cC = mesh.cellCells()[i];
                bool adjToSolid = false;
                for (int j = 0; j < cC.size(); j++)
                {
                    const auto adjCI = cC[j];
                    if (isEqual(marker[adjCI], cSolid))
                    {
                        adjToSolid = true;
                        break;
                    }
                }
                if (!adjToSolid)
                {
                    marker[i] = cFluid;
                    --nGhost;
                    ++nFluid;
                }
            }
        }
    }

    if (nFluid > 0)
    {
        // Remove redundent ghost cells on solid side (through face connection)
        for (int i = 0; i < mesh.nCells(); i++)
        {
            if (isEqual(marker[i], cGhost))
            {
                const auto &cC = mesh.cellCells()[i];
                bool adjToFluid = false;
                for (int j = 0; j < cC.size(); j++)
                {
                    const auto adjCI = cC[j];
                    if (isEqual(marker[adjCI], cFluid))
                    {
                        adjToFluid = true;
                        break;
                    }
                }
                if (!adjToFluid)
                {
                    marker[i] = cSolid;
                    --nGhost;
                    ++nSolid;
                }
            }
        }
    }

    Foam::Info << nFluid << "(fluid) + " << nGhost << "(ghost) + " << nSolid << "(solid) = " << nFluid+nGhost+nSolid << "/" << mesh.nCells() << Foam::endl;
}

void ibInterp_Dirichlet_Linear(const Foam::vector &p_src, Foam::scalar val_src, const std::vector<Foam::vector> &p_adj, const std::vector<Foam::scalar> &val_adj, const Foam::vector &p_dst, Foam::scalar &val_dst)
{
    const int nADJ = p_adj.size();
    double dx, dz;

    Eigen::MatrixXd A(nADJ, 2);
    Eigen::VectorXd b(nADJ);
    Eigen::VectorXd C(2);

    // Assemble adjacent interpolation equations
    for (int i = 0; i < nADJ; i++)
    {
        b[i] = val_adj[i] - val_src;

        dx = p_adj[i].x() - p_src.x();
        dz = p_adj[i].z() - p_src.z();

        A(i, 0) = dx;
        A(i, 1) = dz;
    }

    // Solve least-squares coefficients
    C = A.colPivHouseholderQr().solve(b);

    // Interpolate at the image point
    const Foam::vector p_img = 2.0 * p_src - p_dst;
    dx = p_img.x() - p_src.x();
    dz = p_img.z() - p_src.z();
    const Foam::scalar val_img = val_src + C[0] * dx + C[1] * dz;

    // Interpolate at the ghost point
    val_dst = 2.0 * val_src - val_img;
}

void ibInterp_Neumann_Linear(const Foam::vector &p_src, const Foam::vector &n_src, Foam::scalar snGrad_src, const std::vector<Foam::vector> &p_adj, const std::vector<Foam::scalar> &val_adj, const Foam::vector &p_dst, Foam::scalar &val_dst)
{
    const int nADJ = p_adj.size();
    double dx, dz;

    Eigen::MatrixXd A(nADJ+1, 3);
    Eigen::VectorXd b(nADJ+1);
    Eigen::VectorXd C(3);
    
    // Assemble adjacent interpolation equations
    for (int i = 0; i < nADJ; i++)
    {
        b[i] = val_adj[i];

        dx = p_adj[i].x() - p_src.x();
        dz = p_adj[i].z() - p_src.z();

        A(i, 0) = 1.0;
        A(i, 1) = dx;
        A(i, 2) = dz;
    }

    // Incoporate the boundary condition
    b[nADJ] = snGrad_src;
    A(nADJ, 0) = 0.0;
    A(nADJ, 1) = n_src.x();
    A(nADJ, 2) = n_src.z();

    // Solve least-squares coefficients
    C = A.colPivHouseholderQr().solve(b);

    // Interpolate at the ghost point
    dx = p_dst.x() - p_src.x();
    dz = p_dst.z() - p_src.z();
    val_dst = C[0] + C[1] * dx + C[2] * dz;
}

void ibInterp_zeroGradient_Linear(const Foam::vector &p_src, const Foam::vector &n_src, const std::vector<Foam::vector> &p_adj, const std::vector<Foam::scalar> &val_adj, const Foam::vector &p_dst, Foam::scalar &val_dst)
{
    const int nADJ = p_adj.size();
    double dx, dz;

    Eigen::MatrixXd A(nADJ+1, 3);
    Eigen::VectorXd b(nADJ+1);
    Eigen::VectorXd C(3);
    
    // Assemble adjacent interpolation equations
    for (int i = 0; i < nADJ; i++)
    {
        b[i] = val_adj[i];

        dx = p_adj[i].x() - p_src.x();
        dz = p_adj[i].z() - p_src.z();

        A(i, 0) = 1.0;
        A(i, 1) = dx;
        A(i, 2) = dz;
    }

    // Incoporate the boundary condition
    b[nADJ] = 0.0;
    A(nADJ, 0) = 0.0;
    A(nADJ, 1) = n_src.x();
    A(nADJ, 2) = n_src.z();

    // Solve least-squares coefficients
    C = A.colPivHouseholderQr().solve(b);

    // Interpolate at the image point
    const Foam::vector p_img = 2.0 * p_src - p_dst;
    dx = p_img.x() - p_src.x();
    dz = p_img.z() - p_src.z();
    const Foam::scalar val_img = C[0] + C[1] * dx + C[2] * dz;

    // Interpolate at the ghost point
    val_dst = val_img;
}

void ibInterp_Robin_Linear(const Foam::vector &p_src, const Foam::vector &n_src, Foam::scalar bc_alpha, Foam::scalar bc_beta, Foam::scalar bc_gamma, const std::vector<Foam::vector> &p_adj, const std::vector<Foam::scalar> &val_adj, const Foam::vector &p_dst, Foam::scalar &val_dst)
{
    const int nADJ = p_adj.size();
    double dx, dz;

    Eigen::MatrixXd A(nADJ+1, 3);
    Eigen::VectorXd b(nADJ+1);
    Eigen::VectorXd C(3);

    // Assemble adjacent interpolation equations
    for (int i = 0; i < nADJ; i++)
    {
        b[i] = val_adj[i];

        dx = p_adj[i].x() - p_src.x();
        dz = p_adj[i].z() - p_src.z();

        A(i, 0) = 1.0;
        A(i, 1) = dx;
        A(i, 2) = dz;
    }

    // Incoporate the boundary condition
    b[nADJ] = bc_gamma;
    A(nADJ, 0) = bc_beta;
    A(nADJ, 1) = bc_alpha * n_src.x();
    A(nADJ, 2) = bc_alpha * n_src.z();

    // Solve least-squares coefficients
    C = A.colPivHouseholderQr().solve(b);

    // Interpolate at the image point
    const Foam::vector p_img = 2.0 * p_src - p_dst;
    dx = p_img.x() - p_src.x();
    dz = p_img.z() - p_src.z();
    const double val_img = C[0] + C[1] * dx + C[2] * dz;

    // Interpolate at the ghost point
    double snGrad = C[1] * n_src.x() + C[2] * n_src.z();
    double d = Foam::mag(p_dst - p_img);
    val_dst = val_img - snGrad * d;
}

void ibInterp_Dirichlet_Quadratic(const Foam::vector &p_src, Foam::scalar val_src, const std::vector<Foam::vector> &p_adj, const std::vector<Foam::scalar> &val_adj, const Foam::vector &p_dst, Foam::scalar &val_dst)
{
    const int nADJ = p_adj.size();
    double dx, dz;

    Eigen::MatrixXd A(nADJ, 5);
    Eigen::VectorXd b(nADJ);
    Eigen::VectorXd C(5);

    for (int i = 0; i < nADJ; i++)
    {
        b[i] = val_adj[i] - val_src;

        dx = p_adj[i].x() - p_src.x();
        dz = p_adj[i].z() - p_src.z();

        A(i, 0) = dx;
        A(i, 1) = dz;
        A(i, 2) = dx*dz;
        A(i, 3) = dx*dx;
        A(i, 4) = dz*dz;
    }

    C = A.colPivHouseholderQr().solve(b);

    dx = p_dst.x() - p_src.x();
    dz = p_dst.z() - p_src.z();
    val_dst = val_src + C[0] * dx + C[1] * dz + C[2] * dx*dz + C[3] * dx*dx + C[4] * dz*dz;
}

void normalizeY(double &Y1, double &Y2, double &Y3, double &Y4)
{
    Y1 = std::max(0.0, std::min(Y1, 1.0));
    Y2 = std::max(0.0, std::min(Y2, 1.0));
    Y3 = std::max(0.0, std::min(Y3, 1.0));
    Y4 = std::max(0.0, std::min(Y4, 1.0));

    const double Ysum = Y1+Y2+Y3+Y4;
    Y1 /= Ysum;
    Y2 /= Ysum;
    Y3 /= Ysum;
    Y4 /= Ysum;
}

void calc_gas_property(double T, double &lambda_, double &mu_, double &Cp_)
{
    lambda_ = 1.08e-4 * T + 0.0133;              // Unit: W/m/K
    Cp_ = 0.3 * kcal2J;                          // Unit: J/kg/K
    mu_ = Pr * lambda_ / Cp_;                    // Unit: kg/m/s
}

/**
 * Chemical source terms calculated from the finite-rate global reactions.
 * ----------------------------------------------
 *     Y1 -> Y3                      R1 (Qg1)
 *     beta*Y1 + Y2 -> (beta+1)*Y4   R2 (Qg2)
 *     beta*Y3 + Y2 -> (beta+1)*Y4   R3 (Qg3)
 * ----------------------------------------------
 * @param T Temperature, Unit: K
 * @param p Pressure, Unit: bar
 * @param Y1 Mass fraction of the 1st symbolic species
 * @param Y2 Mass fraction of the 2nd symbolic species
 * @param Y3 Mass fraction of the 3rd symbolic species
 * @param wdot1 Production rate of the 1st symbolic species, Unit: kg/m^3/s
 * @param wdot2 Production rate of the 2nd symbolic species, Unit: kg/m^3/s
 * @param wdot3 Production rate of the 3rd symbolic species, Unit: kg/m^3/s
 * @param wdot4 Production rate of the 4th symbolic species, Unit: kg/m^3/s
 * @param Qdot Heat source due to chemical reactions, Unit: W/m^3
 */
void calc_chem_src(double T, double p, double Y1, double Y2, double Y3, double &wdot1, double &wdot2, double &wdot3, double &wdot4, double &Qdot)
{
    // Reaction rates, Unit: kg/m^3/s
    const double R1 = 4.11e4 * std::pow(p, 2.06) * Y1 * std::exp(-3000.0/T);
    const double R2 = 2.35e7 * std::pow(p, 2.06) * std::pow(Y1, 3.3) * std::pow(Y2, 0.4) * std::exp(-8500.0/T);
    const double R3 = 9.5e4 * std::pow(p, 1.6) * Y2 * Y3 * std::exp(-8500.0/T);

    // Source terms
    const double beta = 6.89;
    wdot1 = -R1 - beta * R2;
    wdot2 = -R1 - R2;
    wdot3 = R1 - beta * R3;
    wdot4 = (beta+1) * (R2 + R3);
    Qdot = (410.0 * R1 + 7403.0 * R2 + 4396.0 * R3) * kcal2J;
}

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"

    /* Load meshes */
    #include "createMeshes.H"

    /* Create variables */
    #include "createFields.H"

    /* Initialize */
    {
        // Density
        rho = p / Rg / T;
        rho.correctBoundaryConditions();

        // Signed-distance
        for (int i = 0; i < mesh_gas.nCells(); i++)
        {
            phi[i] = mesh_gas.C()[i].z() - plane_z;
        }
        for (int i = 0; i < pointMesh_gas.size(); i++)
        {
            phi_pnt[i] += mesh_gas.points()[i].z() - plane_z;
        }

        // Properties
        for (int i = 0; i < mesh_gas.nCells(); i++)
        {
            calc_gas_property(T[i], lambda[i], mu[i], Cp[i]);
        }
        mu.correctBoundaryConditions();
        Cp.correctBoundaryConditions();
        lambda.correctBoundaryConditions();
    }

    /* createPhi.H */
    USn = Foam::fvc::interpolate(U) & mesh_gas.Sf();
    setBdryVal(mesh_gas, U, USn);

    /* compressibleCreatePhi.H */
    rhoUSn = Foam::fvc::interpolate(rho*U) & mesh_gas.Sf();
    setBdryVal(mesh_gas, rho, U, rhoUSn);

    /* Record adjacent cells of each cell (degree 2) */
    Foam::List<adj2Stencil> adjInfo(mesh_gas.nCells());
    for (int i = 0; i < mesh_gas.nCells(); i++)
    {
        Foam::labelHashSet adj;

        for (int j = 0; j < mesh_gas.cellPoints()[i].size(); j++)
        {
            const auto pI = mesh_gas.cellPoints()[i][j];
            for (int k = 0; k < mesh_gas.pointCells()[pI].size(); k++)
            {
                const auto cI = mesh_gas.pointCells()[pI][k];
                if (cI != i)
                    adj.insert(cI);
            }
        }

        for (int j = 0; j < mesh_gas.cellCells()[i].size(); j++)
        {
            const auto cI = mesh_gas.cellCells()[i][j];
            for (int k = 0; k < mesh_gas.cellCells()[cI].size(); k++)
            {
                const auto cI2 = mesh_gas.cellCells()[cI][k];
                if (cI2 != i)
                    adj.insert(cI2);
            }
        }

        adjInfo[i].myC = i;
        adjInfo[i].adjC = adj.toc();
    }

    while(runTime.loop())
    {
        const Foam::dimensionedScalar dt(Foam::dimTime, runTime.deltaTValue());
        Foam::Info << "\nn=" << runTime.timeIndex() << ", t=" << std::stod(runTime.timeName(), nullptr)*s2ms << "ms, dt=" << dt.value()*s2ns << "ns" << Foam::endl;
        runTime.write();

        /* Classify mesh cells & interpolate on newly exposed fresh cells */
        if (runTime.timeIndex() == 1)
        {
            identifyIBCell(mesh_gas, phi, phi_pnt, cIbMarker);
            for (int i = 0; i < mesh_gas.nCells(); i++)
                cIbMask[i] = cIbMarker[i] > 0 ? 1.0 : 0.0;
        }
        else
        {
            const Foam::volScalarField cIbMarker_old(cIbMarker);

            identifyIBCell(mesh_gas, phi, phi_pnt, cIbMarker);
            for (int i = 0; i < mesh_gas.nCells(); i++)
                cIbMask[i] = cIbMarker[i] > 0 ? 1.0 : 0.0;

            // Clear cell fresh flag
            for (int i = 0; i < mesh_gas.nCells(); i++)
                cIbFresh[i] = 0.0;

            // Record and mark fresh cells
            Foam::scalarList fresh;
            for (int i = 0; i < mesh_gas.nCells(); i++)
            {
                if(!isEqual(cIbMarker_old[i], cFluid) && isEqual(cIbMarker[i], cFluid))
                {
                    fresh.append(i);
                    cIbFresh[i] = cFresh;
                }
            }
            if (!fresh.empty())
                Foam::Info << fresh.size() << " fresh cells to be handled!" << Foam::endl;

            // Interpolate on fresh cells
            for (int i = 0; i < fresh.size(); i++)
            {
                const auto cI = fresh[i];

                Foam::Info << "Cell[" << cI << "]: rho=" << rho[cI] << ", U=" << U[cI] << ", p=" << p[cI] << ", T=" << T[cI] << ", Y1=" << Y1[cI] << ", Y2=" << Y2[cI] << ", Y3=" << Y3[cI] << ", Y4=" << Y4[cI] << Foam::endl;

                // Collect adjacent data
                Foam::labelList idx;
                for (int j = 0; j < adjInfo[cI].adjC.size(); j++)
                {
                    const auto adjCI = adjInfo[cI].adjC[j];
                    if (isEqual(cIbMarker[adjCI], cFluid) && !isEqual(cIbFresh[adjCI], cFresh))
                        idx.append(adjCI);
                }
                if (idx.size() < 2)
                {
                    Foam::Info << "Too less adjacent cells! (" << idx.size() << ") of cell[" << cI << "]" << Foam::endl;

                    // Foam::Info << ">>>>>>>>>>Debugging>>>>>>>>>>" << Foam::endl;
                    // for (int j = 0; j < adjInfo[cI].adjC.size(); j++)
                    // {
                    //     const auto adjCI = adjInfo[cI].adjC[j];
                    //     Foam::Info << "Adjacent cell[" << adjCI << "] (" << j << "): cIbMarker=" << cIbMarker[adjCI] << ", cIbFresh=" << cIbFresh[adjCI] << Foam::endl;
                    // }
                    // Foam::Info << "<<<<<<<<<<END<<<<<<<<<<" << Foam::endl;

                    return -2;
                }

                Foam::Info << "Cell[" << cI << "]: " << idx.size() << " adjacent cells for interpolation" << Foam::endl;

                std::vector<Foam::scalar> val_ux(idx.size()), val_uz(idx.size());
                std::vector<Foam::scalar> val_p(idx.size());
                std::vector<Foam::scalar> val_T(idx.size());
                std::vector<Foam::scalar> val_Y1(idx.size()), val_Y2(idx.size()), val_Y3(idx.size()), val_Y4(idx.size());
                std::vector<Foam::vector> pos(idx.size());
                for (int j = 0; j < idx.size(); j++)
                {
                    const auto cI = idx[j];
                    val_ux[j] = U[cI].x();
                    val_uz[j] = U[cI].z();
                    val_p[j] = p[cI];
                    val_T[j] = T[cI];
                    val_Y1[j] = Y1[cI];
                    val_Y2[j] = Y2[cI];
                    val_Y3[j] = Y3[cI];
                    val_Y4[j] = Y4[cI];
                    pos[j] = mesh_gas.C()[cI];
                }

                // Boundary-Intersection position
                const Foam::vector p_BI(mesh_gas.C()[cI].x(), 0.0, mesh_gas.C()[cI].z()-phi[cI]);
                const Foam::vector n_BI(0.0, 0.0, 1.0);

                // Velocity
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_ux, mesh_gas.C()[cI], U[cI].x());
                ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_uz, mesh_gas.C()[cI], U[cI].z());
                U[cI].y() = 0.0;

                // Pressure
                ibInterp_zeroGradient_Linear(p_BI, n_BI, pos, val_p, mesh_gas.C()[cI], p[cI]);

                // Temperature
                ibInterp_Dirichlet_Linear(p_BI, plane_T, pos, val_T, mesh_gas.C()[cI], T[cI]);

                // Species
                if (mesh_gas.C()[cI].x() < 0.3*L || mesh_gas.C()[cI].x() > 0.7*L)
                {
                    ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_Y1, mesh_gas.C()[cI], Y1[cI]);
                    ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y2, mesh_gas.C()[cI], Y2[cI]);
                }
                else
                {
                    ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y1, mesh_gas.C()[cI], Y1[cI]);
                    ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_Y2, mesh_gas.C()[cI], Y2[cI]);
                }
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y3, mesh_gas.C()[cI], Y3[cI]);
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y4, mesh_gas.C()[cI], Y4[cI]);

                // Density
                rho[cI] = p[cI] / Rg.value() / T[cI];

                Foam::Info << "Cell[" << cI << "]: rho=" << rho[cI] << ", U=" << U[cI] << ", p=" << p[cI] << ", T=" << T[cI] << ", Y1=" << Y1[cI] << ", Y2=" << Y2[cI] << ", Y3=" << Y3[cI] << ", Y4=" << Y4[cI] << Foam::endl;
            }

            rhoUSn = Foam::fvc::interpolate(rho*U) & mesh_gas.Sf();
            setBdryVal(mesh_gas, rho, U, rhoUSn);
        }

        /* Interpolation on IB cells */
        for (int i = 0; i < mesh_gas.nCells(); i++)
        {
            if (isEqual(cIbMarker[i], cGhost))
            {
                // Extract data on neighboring fluid cells
                std::vector<Foam::label> idx;
                for (int j = 0; j < adjInfo[i].adjC.size(); j++)
                {
                    const auto cI = adjInfo[i].adjC[j];
                    if (isEqual(cIbMarker[cI], cFluid))
                        idx.push_back(cI);
                }
                std::vector<Foam::scalar> val_ux(idx.size()), val_uz(idx.size());
                std::vector<Foam::scalar> val_p(idx.size());
                std::vector<Foam::scalar> val_T(idx.size());
                std::vector<Foam::scalar> val_Y1(idx.size()), val_Y2(idx.size()), val_Y3(idx.size()), val_Y4(idx.size());
                std::vector<Foam::vector> pos(idx.size());
                for (int j = 0; j < idx.size(); j++)
                {
                    const auto cI = idx[j];
                    val_ux[j] = U[cI].x();
                    val_uz[j] = U[cI].z();
                    val_p[j] = p[cI];
                    val_T[j] = T[cI];
                    val_Y1[j] = Y1[cI];
                    val_Y2[j] = Y2[cI];
                    val_Y3[j] = Y3[cI];
                    val_Y4[j] = Y4[cI];
                    pos[j] = mesh_gas.C()[cI];
                }

                // Boundary-Intersection position
                const Foam::vector p_BI(mesh_gas.C()[i].x(), 0.0, mesh_gas.C()[i].z()-phi[i]);
                const Foam::vector n_BI(0.0, 0.0, 1.0);

                // Velocity
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_ux, mesh_gas.C()[i], U[i].x());
                ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_uz, mesh_gas.C()[i], U[i].z());
                U[i].y() = 0.0;

                // Pressure
                ibInterp_zeroGradient_Linear(p_BI, n_BI, pos, val_p, mesh_gas.C()[i], p[i]);

                // Temperature
                ibInterp_Dirichlet_Linear(p_BI, plane_T, pos, val_T, mesh_gas.C()[i], T[i]);

                // Species
                if (mesh_gas.C()[i].x() < 0.3*L || mesh_gas.C()[i].x() > 0.7*L)
                {
                    ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_Y1, mesh_gas.C()[i], Y1[i]);
                    ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y2, mesh_gas.C()[i], Y2[i]);
                }
                else
                {
                    ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y1, mesh_gas.C()[i], Y1[i]);
                    ibInterp_Dirichlet_Linear(p_BI, 1.0, pos, val_Y2, mesh_gas.C()[i], Y2[i]);
                }
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y3, mesh_gas.C()[i], Y3[i]);
                ibInterp_Dirichlet_Linear(p_BI, 0.0, pos, val_Y4, mesh_gas.C()[i], Y4[i]);
            }
        }

        /* Update reaction rates */
        for (int i = 0; i < mesh_gas.nCells(); i++)
        {
            if (isEqual(cIbMarker[i], cFluid))
            {
                if (runTime.timeIndex() == 17701)
                    Foam::Info << "cell[" << i << "]: rho=" << rho[i] << ", U=" << U[i] << ", p=" << p[i]/one_bar << ", T=" << T[i] << ", Y1=" << Y1[i] << ", Y2=" << Y2[i] << ", Y3=" << Y3[i] << ", Y4=" << Y4[i] << Foam::endl;

                calc_chem_src(T[i], p[i]/one_bar, Y1[i], Y2[i], Y3[i], wdot1[i], wdot2[i], wdot3[i], wdot4[i], Qdot[i]);
            }
            else
            {
                wdot1[i] = wdot2[i] = wdot3[i] = wdot4[i] = 0.0;
                Qdot[i] = 0.0;
            }
        }

        /* Update gas-phase by semi-implicit iteration */
        bool converged = false;
        int m = 0;
        while(++m <= 5)
        {
            Foam::Info << "m=" << m << Foam::endl;

            /* Predictor */
            {
                /* Provisional density */
                {
                    Foam::fvScalarMatrix rhoEqn
                    (
                        Foam::fvm::ddt(rho) + Foam::fvc::div(rhoUSn)
                    );
                    rhoEqn.solve();
                    rho_star = rho;
                }

                /* Temperature */
                {
                    Foam::fvScalarMatrix TEqn
                    (
                        Foam::fvm::ddt(rho, T) + Foam::fvc::div(rhoUSn, T) == (Foam::fvc::laplacian(lambda, T) + Qdot) / Cp
                    );

                    // For solid cells, set temperature to constant;
                    // For ghost cells, set to interpolated value.
                    {
                        Foam::scalarList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(T0);
                                else
                                    val.append(std::max(T0, T[i]));
                            }
                        }
                        TEqn.setValues(idx, val);
                    }

                    TEqn.solve();
                }

                /* Species */
                {
                    Foam::fvScalarMatrix Y1Eqn
                    (
                        Foam::fvm::ddt(rho, Y1) + Foam::fvc::div(rhoUSn, Y1) == Foam::fvc::laplacian(lambda/Cp, Y1) + wdot1
                    );

                    // For solid cells, set species mass fraction to zero;
                    // For ghost cells, set to interpolated value.
                    {
                        Foam::scalarList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(0.0);
                                else
                                    val.append(Y1[i]);
                            }
                        }
                        Y1Eqn.setValues(idx, val);
                    }

                    // Solve
                    Y1Eqn.solve();
                }
                {
                    Foam::fvScalarMatrix Y2Eqn
                    (
                        Foam::fvm::ddt(rho, Y2) + Foam::fvc::div(rhoUSn, Y2) == Foam::fvc::laplacian(lambda/Cp, Y2) + wdot2
                    );

                    // For solid cells, set species mass fraction to zero;
                    // For ghost cells, set to interpolated value.
                    {
                        Foam::scalarList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(0.0);
                                else
                                    val.append(Y2[i]);
                            }
                        }
                        Y2Eqn.setValues(idx, val);
                    }

                    // Solve
                    Y2Eqn.solve();
                }
                {
                    Foam::fvScalarMatrix Y3Eqn
                    (
                        Foam::fvm::ddt(rho, Y3) + Foam::fvc::div(rhoUSn, Y3) == Foam::fvc::laplacian(lambda/Cp, Y3) + wdot3
                    );

                    // For solid cells, set species mass fraction to zero;
                    // For ghost cells, set to interpolated value.
                    {
                        Foam::scalarList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(0.0);
                                else
                                    val.append(Y3[i]);
                            }
                        }
                        Y3Eqn.setValues(idx, val);
                    }

                    // Solve
                    Y3Eqn.solve();
                }
                {
                    Foam::fvScalarMatrix Y4Eqn
                    (
                        Foam::fvm::ddt(rho, Y4) + Foam::fvc::div(rhoUSn, Y4) == Foam::fvc::laplacian(lambda/Cp, Y4) + wdot4
                    );

                    // For solid cells, set species mass fraction to zero;
                    // For ghost cells, set to interpolated value.
                    {
                        Foam::scalarList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(0.0);
                                else
                                    val.append(Y4[i]);
                            }
                        }
                        Y4Eqn.setValues(idx, val);
                    }

                    // Solve
                    Y4Eqn.solve();
                }

                /* Ensure the sum of Y being unity */
                for (int i=0; i < mesh_gas.nCells(); i++)
                {
                    if (cIbMarker[i] > cSolid)
                        normalizeY(Y1[i], Y2[i], Y3[i], Y4[i]);
                }
                Y1.correctBoundaryConditions();
                Y2.correctBoundaryConditions();
                Y3.correctBoundaryConditions();
                Y4.correctBoundaryConditions();

                /* Provisional momentum */
                {
                    grad_p = Foam::fvc::grad(p);
                    Foam::fvVectorMatrix UEqn
                    (
                        Foam::fvm::ddt(rho, U) + Foam::fvc::div(rhoUSn, U) == -grad_p + Foam::fvc::laplacian(mu, U)
                    );

                    // For solid cells, set velocity to zero;
                    // For ghost cells, set velocity to interpolated value.
                    {
                        Foam::vectorList val;
                        Foam::labelList idx;

                        for (int i=0; i < mesh_gas.nCells(); i++)
                        {
                            if (cIbMarker[i] < cFluid)
                            {
                                idx.append(i);
                                if (isEqual(cIbMarker[i], cSolid))
                                    val.append(Foam::vector(0.0, 0.0, 0.0));
                                else
                                    val.append(U[i]);
                            }
                        }
                        UEqn.setValues(idx, val);
                    }

                    UEqn.solve();
                }

                /* Provisional mass flux */
                rhoUSn = Foam::fvc::interpolate(rho_star*U) & mesh_gas.Sf();
                Foam::surfaceScalarField grad_p_f_compact(Foam::fvc::snGrad(p) * mesh_gas.magSf());
                Foam::surfaceScalarField grad_p_f_mean(Foam::fvc::interpolate(grad_p) & mesh_gas.Sf());
                rhoUSn -= dt * (grad_p_f_compact - grad_p_f_mean);
                setBdryVal(mesh_gas, rho, U, rhoUSn);

                /* Density (thermal) */
                rho = p / Rg / T;
                rho.correctBoundaryConditions();
            }

            /* Corrector */
            {
                /* Pressure-correction Helmholtz equation */
                Foam::fvScalarMatrix dpEqn
                (
                    Foam::fvm::Sp(1.0/(dt * Rg * T), dp) - dt * Foam::fvm::laplacian(dp) == -(Foam::fvc::div(rhoUSn) + Foam::fvc::ddt(rho))
                );

                // For both solid and ghost cells, set the incremental pressure-correction to zero.
                {
                    Foam::scalarList val;
                    Foam::labelList idx;

                    for (int i=0; i < mesh_gas.nCells(); i++)
                    {
                        if (cIbMarker[i] < cFluid)
                        {
                            idx.append(i);
                            val.append(0.0);
                        }
                    }
                    dpEqn.setValues(idx, val);
                }

                dpEqn.solve();

                /* Update */
                p += dp;
                p.correctBoundaryConditions();

                U = (rho_star*U - dt * Foam::fvc::grad(dp)) / rho;
                U.correctBoundaryConditions();

                rhoUSn -= dt * Foam::fvc::snGrad(dp) * mesh_gas.magSf();
                setBdryVal(mesh_gas, rho, U, rhoUSn);

                USn = rhoUSn / Foam::fvc::interpolate(rho);
                setBdryVal(mesh_gas, U, USn);
            }

            /* Check convergence */
            {
                Foam::scalar eps_inf, eps_1, eps_2;

                diagnose(mesh_gas, meshInfo_gas, dp, cIbMask, eps_1, eps_2, eps_inf);
                Foam::Info << "||dp||: " << eps_inf << "(Inf), " << eps_1 << "(1), " << eps_2 << "(2)" << Foam::endl;
                const bool criteria_dp = eps_inf < 1e-5 * p0 && eps_1 < 1e-6 * p0 && eps_2 < 1e-5 * p0;

                diagnose(mesh_gas, meshInfo_gas, rho-rho_star, cIbMask, eps_1, eps_2, eps_inf);
                Foam::Info << "||rho(m)-rho*||: " << eps_inf << "(Inf), " << eps_1 << "(1), " << eps_2 << "(2)" << Foam::endl;
                const bool criteria_drho = eps_inf < 1e-3 || eps_1 < 1e-6 || eps_2 < 1e-5;

                converged = criteria_dp && criteria_drho;
                if (converged && m > 1)
                    break;
            }
        }
        if (!converged)
            Foam::Info << "\nGas-phase failed to converged after " << m-1 << " semi-implicit iterations!" << Foam::endl;

        /* Update gas-phase properties */
        for (int i = 0; i < mesh_gas.nCells(); i++)
        {
            calc_gas_property(T[i], lambda[i], mu[i], Cp[i]);
        }
        mu.correctBoundaryConditions();
        Cp.correctBoundaryConditions();
        lambda.correctBoundaryConditions();

        /* Check range of gas-phase variables */
        {
            Foam::scalar eps_inf, eps_1, eps_2;
            Foam::scalar vMin, vMax;

            // Density
            diagnose(mesh_gas, rho, cIbMask, vMin, vMax);
            Foam::Info << "\nrho: " << vMin << " ~ " << vMax << Foam::endl;

            // Velocity magnitude
            diagnose(mesh_gas, Foam::mag(U), cIbMask, vMin, vMax);
            Foam::Info << "|U|: " << vMin << " ~ " << vMax << Foam::endl;

            // Pressure fluctuation
            diagnose(mesh_gas, p, cIbMask, vMin, vMax);
            Foam::Info << "p-p0: " << vMin-p0 << " ~ " << vMax-p0 << Foam::endl;

            // Temperature
            diagnose(mesh_gas, T, cIbMask, vMin, vMax);
            Foam::Info << "T: " << vMin << " ~ " << vMax << Foam::endl;
            // if (vMax > 3500.0 || vMin < 100.0)
            // {
            //     Foam::Info << "Invalid temperature detected!" << Foam::endl;
            //     return -3;
            // }

            // Continuity
            diagnose(mesh_gas, meshInfo_gas, Foam::fvc::div(rhoUSn), cIbMask, eps_1, eps_2, eps_inf);
            Foam::Info << "||div(rhoU)||: " << eps_inf << "(Inf), " << eps_1 << "(1), " << eps_2 << "(2)" << Foam::endl;
        }

        /* Update the position of gas-solid interface to (n+1) */
        {
            for (int i = 0; i < mesh_gas.nCells(); i++)
                phi[i] += plane_rb * dt.value();
            phi.correctBoundaryConditions();

            for (int i = 0; i < pointMesh_gas.size(); i++)
                phi_pnt[i] += plane_rb * dt.value();
            phi_pnt.correctBoundaryConditions();
        }
    }

    return 0;
}
